Initial data for black hole— neutron star binaries: a flexible, high-accuracy spectral 

method. 
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We present a new numerical scheme to solve the initial value problem for black hole-neutron 
star binaries. This method takes advantage of the flexibility and fast convergence of a multidomain 
spectral representation of the initial data to construct high-accuracy solutions at a relatively low 
computational cost. We provide convergence tests of the method for both isolated neutron stars 
and irrotational binaries. In the second case, we show that we can resolve the small inconsistencies 
that are part of the quasiequilibrium formulation, and that these inconsistencies are significantly 
smaller than observed in previous works. The possibility of generating a wide variety of initial data 
is also demonstrated through two new configurations inspired by results from binary black holes. 
First, we show that choosing a modified Kerr-Schild conformal metric instead of a flat conformal 
metric allows for the construction of quasiequilibrium binaries with a spinning black hole. Second, 
we construct binaries in low-eccentricity orbits, which are a better approximation to astrophysical 
binaries than quasiequilibrium systems. 

PACS numbers: 04.25.dk,04.40.Dg,04.30.Db,04.20.Ex,95.30.Sf 



I. INTRODUCTION 

Over the last few years, the prospect of gravitational 
wave detection by ground based experiments such as 
LIGO [l| and VIRGO Q has encouraged rapid de- 
velopments in the field of numerical relativity. Most of 
that effort was aimed at the evolution of compact bi- 
naries, sources of waves potentially observable by those 
detectors. Binary neutron stars were the first to be suc- 
cessfully evolved in a fully relativistic framework, and 
have been studied regularly over the last eight years 
S i, H, i, S S, il- Evolutions of binaryblack holes 
(BBH) followed a few years later [10, [HI, [l^l ■, ^tnd con- 
tinue to be an extraordinarily active area of research (see 
[l3t and references therein). 

The third type of compact binary, black hole-neutron 
star (BH-NS) binaries, has not been as widely studied 
yet. The evolution of the black hole singularity and the 
presence of matter combine the difficulties of evolving 
both binary black holes and binary neutron stars. And 
the system has its own specific challenges, notably the 
accretion of the neutron star matter onto the black hole. 
Such binaries are, however, worth studying not only for 
their interest as gravitational wave sources, but also as 
potential sources of gamma-ray bursts [1J|. The first 
evolutions of such systems were announced very recently 
jlSl . Il6[, and such evolutions will allow more extensive 
study of their wave emission, merger, accretion disk for- 
mation, and so on. 

The choice of a suitable initial configuration for bi- 
nary evolutions has been a long-standing problem. Not 
only do the Einstein equations include constraints on the 
initial data, but also choosing a starting point that rep- 
resents a realistic astrophysical situation is not trivial. 
Because of their computational cost, numerical simula- 
tions of compact binaries usually start just a few orbits 



away from merger. The two objects are close enough 
that the nonlinearity of the Einstein equations is impor- 
tant. In that regime, there is no known way of prescrib- 
ing the exact state of the system. The most common 
assumption is that the binary has had time to settle into 
a quasiequilibrium state, the system being approximately 
time-independent in the corotating frame. Furthermore, 
as the viscous forces within the star are expected to be 
small, we do not expect much change in the spin of the 
star as the orbital radius decreases. For an initially non- 
spinning neutron star, this would lead to an irrotational 
velocity profile, another standard assumption. Because 
of gravitational wave emission, there is no exact equilib- 
rium state, however. Accordingly, these conditions can- 
not be perfectly satisfied, a problem we will discuss in 
more detail later on. 

Previous results on initial data for BH-NS evolutions 
include the early work of Taniguchi et al. [131 ^^^ ?>0Y>- 
uerta et al. [T^], as well as more recent initial config- 
urations generated by Taniguchi et al. [l^, [20, [2l| and 
Grandclement [2^. Both Taniguchi and Grandclement 
use codes based on the lorene package [231, and their 
most recent publications are similar in accuracy, compu- 
tational cost, and numerical results. 

In this paper, we present an alternative numerical 
scheme for the solution of this problem. Our code is 
based on the spectral elliptic solver (spells) developed 
by the Cornell-Caltech collaboration [24|, and originally 
used by Pfeiffer [2^, [20] for the study of binary black 
holes (BBH) initial data. For our numerical tests, the 
mathematical formulation of the problem will be very 
similar to [21[ and [22| , allowing easy evaluation of the 
performance of our code. 

Our motivation for using spells is the remarkable flex- 
ibility of its multidomain spectral methods. This allows 
us to efficiently adapt the configuration of our numeri- 



cal grid to the geometry of the system and yields high- 
precision results at a very reasonable computational cost. 
As we will see in Sec. lU elliptic equations form the core 
of the initial data problem. Using lorene, each of those 
equations has to be approximated by two Poisson equa- 
tions, with coupled source terms. These two sets of equa- 
tions are then solved through an iterative method. The 
variables are fields with an approximate spherical sym- 
metry around one of the compact objects. However, as 
the source terms for the BH fields include terms centered 
around the NS, obtaining high-precision initial data re- 
quires a large angular resolution. 

With SPELLS, by contrast we do not have to limit our- 
selves to spheres around the compact objects. We can 
instead choose among a wide variety of subdomain ge- 
ometries and coordinate mappings. As the basis func- 
tions of our spectral expansion are more adapted to the 
geometry of the solution, a significantly smaller number 
of collocation points are necessary to reach a given accu- 
racy. 

In Sec. lIVBi we will see that the main sources of error 
in our initial data are the approximations introduced by 
the quasiequilibrium formulation. Using spells, we can 
rapidly solve the initial data problem for a large variety of 
configurations to a precision allowing us to resolve these 
errors. We will show that they appear to be significantly 
lower than quoted in [2l|, I24I. For the closest binaries, 
when the distortion of the star limits the precision of any 
spectral method, such precision is no longer possible — 
at least using our current numerical grid. But our error 
remains reasonable, reaching the level of the deviations 
from equilibrium mentioned in l21[ for the most extreme 
cases. 

In addition to the high-precision initial data our re- 
sults provide for evolutions of BH-NS binaries, they 
should also make it possible to explore the limits of the 
quasiequilibrium formalism. Such studies are already 
possible for BBH binaries, as shown in [23|. On the 
BH horizon, the deviations from equilibrium computed 
in 



27[ are similar to our own results. 



Using SPELLS, we are also able to study initial data for 
a spinning BH by abandoning the assumption of confor- 
mal flatness. Earlier results showed that initial configura- 
tions built using a Kerr-Schild conformal metric were sig- 
nificantly inferior to their conformally flat counterparts 
[1^, l20l- Here, adapting a method developed by Lovelace 
for BBH [23, we show that a modified Kerr-Schild met- 
ric can lead to high-precision initial data. In Sec. IIVDI 
we present our results for spinning and nonspinning black 
holes using this modified Kerr-Schild conformal metric. 

We review the formulation of the initial value problem 
in Sec.|TTl and present in more detail our numerical meth- 
ods in Sec. IIIII Then, in Sec. IIVI we discuss some tests 
of our code, including isolated stars and binaries that are 
directly comparable to previous results. Through conver- 
gence tests, we obtain a good estimate of the amplitude of 
constraint violations and of our error in global quantities 
such as the ADM (Arnowitt-Dcscr-Misncr) energy and 



linear and angular momentum. Such convergence tests 
for fully consistent initial data in the presence of matter 
have, to our knowledge, only been published previously in 
the case of NS-NS binaries (see for example [29| , specif- 
ically Figs. 4 to 7), and up to relative precisions slightly 
better than 10^^. Our estimates will confirm that we are 
able to resolve deviations from quasiequilibrium except 
for strongly distorted stars. 

Finally, adapting a method developed by Pfeiffer et al. 
[30| for BBH binaries, we demonstrate the possibility of 
reducing the eccentricity of the system, leading to ini- 
tial configurations more realistic than quasiequilibrium 
orbits. 



II. THE INITIAL DATA PROBLEM 

The construction of initial data on a spatial slice con- 
taining matter typically involves two types of conditions. 
First, from the Einstein equations we know that any ini- 
tial data will have to satisfy the Hamiltonian and momen- 
tum constraints, which we will write as a set of elliptic 
equations. Second, we want the resulting configuration 
to represent a physically reasonable situation. The mass 
of each object, its spin, their initial separation, and the 
cUipticity of the orbit are all parameters we want to con- 
trol, and the initial state and physical properties of the 
fluid have to be carefully chosen. In this section, we will 
describe the different equations used to enforce those con- 
ditions, and their formulation in our numerical solver. 



A. Constraints 

We impose the constraints on our initial spatial slice 
by solving the extended conformal thin sandwich (XCTS) 
system, a set of 5 elliptic equations based on the confor- 
mal thin sandwich decomposition proposed by York I3ll . 
Here, we start from the formulation used by Pfeiffer [2y| 
for BBH binaries, adding the matter contribution as fixed 
source terms in the XCTS equations. 

The metric tensor is written in its 3-1-1 form: 

= -a^dt^ + -fij {dx' + 13' dt) {dx^ + [3^dt) , (1) 

where a is the lapse, /3' the shift, and 7^ the 3-nietric 
induced on a spatial slice at constant t. The normal n to 
such a slice and the tangent to the coordinate line t are 
then related by 



t^^an^+ZJA", 



(2) 



We treat the matter as a perfect fluid and write the stress- 
energy tensor as 



?)<<. = (p + P) u^j_u^ + Pg^u, 



(3) 



where p is the fluid energy density, P its pressure, and transformations according to the scheme ^ 
u^ its 4-velocity. In practice, we will use projections of 
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where U^ is the fluid 3-velocity in the inertial frame, de- 
flncd in terms of the 4-velocity u, the normal n to the 
spatial slice studied, and the Lorentz factor 7„ as 



(11) 
(12) 
(13) 
(14) 
(15) 
(16) 



Denoting the time derivative of the conformal spatial 
metric by Uij = dtgij, Eqs. ^ and pH)) link A'-' and the 
shift by 
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(17) 



u = 7„(n + U). 



(7) 



If the system is close to equilibrium, it is convenient to 
choose the coordinate system so that dt is an approxi- 
mate Killing vector. We will thus try to solve the system 
in coordinates comoving with the binary. In such a co- 
ordinate system, the shift increases in magnitude with 
the distance from the center of rotation and diverges at 
spatial infinity. This is a difficulty for numerical solvers. 
Furthermore, to control the eccentricity of the binary, we 
choose to give the system an initial radial velocity of the 
form V = dor. This also leads to a diverging term in the 
shift at large distances. 



We thus further decompose the shift vector as 



/3 = /3o + r2xr-|- dor, 



(8) 



where the conformal longitudinal operator L is 






(18) 



The XCTS formulation of the constraints is then a set 
of 5 coupled elliptic equations, with the conformal factor 
0, the densitized lapse a<i> = 50^, and the shift /3 (or, in 
practice, the inertial shift /3o) as variables: 
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(19) 
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+2'K(\r^E = 



where /3o is the shift in the inertial frame and il the 
orbital angular velocity of the system. In practice, we 
solve for /3o instead of /3, as /3o conveniently vanishes at 
spatial infinity. We turn now to the extrinsic curvature, 
defined as 



^ fii^ — 2^n9iii' ^ 



(9) 



v^ (a(f) 



{a(t>') 
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h^K^ 



'A,,^^ (21) 



2TT(t)-^[E + 2S 



-(f {dtK - p^dkK) . 



Here, S, S, and J' determine the matter content of the 
slice, and we are free to choose 7ij, My, K, and dtK. 



where £„ is the Lie derivative along the normal n. In 
the conformal thin sandwich formalism, K^j^u is divided 
into its trace K and trace- free part A^^: 



K'^ = A^' + W'K. 



(10) 



The decomposition is completed by the use of conformal 



^ A conformal transformation of the matter quantities E, 3 and 
J' is necessary for the Hamiltonian constraint to have a unique 
solution 1321 . But different choices for the ratio between confor- 
mal and physical quantities are valid. Our choice of (j>^, which 
differs from [2611 . guarantees that volume integrals of the mat- 
ter terms for fixed E, S and J* are independent of the conformal 
factor (f). Indeed, the physical volume element on the spatial slice 
is dV = (fi^^/^d^x, where 7 is the determinant of the conformal 
metric, and thus f EdV = J Ed^x. The full XCTS system is 
known to have non-unique solutions for vacuum |33l . |34| ; this 
may carry over to space-times with matter, but we have not ob- 
served non-uniqueness in the course of the present work. 



Eqs. ((19)) and ([20)) are the momentum and Hamiltonian 
constraints, while Eq. (|^T)) can be derived from the evo- 
lution equation for K^^ . (For more details on the XCTS 
system, and its derivation, see [35[-) 

For quasiequilibrium initial conditions, a natural 
choice for the free variables is to set the time derivatives 
to zero. The choice of Tjw and K is, however, less obvious. 
Taniguchi et al. [l^, [20| showed that a conformally flat 



Sij) with maximal slicing (K ~ 0) gives 



metric {jij 

good results — better than using a Kerr-Schild back- 
ground at least. For the tests in this paper, we will make 
the same choice. In Sec. lIVDi however, we will show 
that different choices lead to acceptable initial data, and 
make it possible to construct spinning BHs. 



B. Hydrostatic equilibrium 

The initial state of the matter within the neutron star 
is, in general, unknown. However, we can make some rea- 
sonable approximations. First, we will require the fluid to 
be in a state of hydrostatic equilibrium in the comoving 
frame. Following the method described by Gourgoulhon 
et al. [231, we use the first integral of the Euler equation. 
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70 



constant. 



(22) 



where h is the fluid enthalpy and we define the Lorentz 
factors 
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As before, V is the fluid 3-velocity in the inertial frame, 
while C/q is the 3-velocity of a comoving observer. For a 
corotating binary, we simply have [/' = [/q, while for an 
irrotational configuration, there should exist a velocity 
potential ^ [23] such that 



U' 






d,^. 



The equation of continuity is then 



Po 
' h 



0, 



(27) 



(28) 



where po is the baryon density. This is an elliptic equa- 
tion in "f, which we can rewrite more explicitly in our 
variables as 



hP'cl)* 






- fd. 



In 



h 



dk^ (29) 



a 



dan + hK-/n4> 



f^d^^djpo 



/»7«/?V' 



diPo- 



For a star in a binary, the main contribution to the po- 
tential ^ comes from the movement of the star along its 
orbit. It is thus convenient to decompose ^ as proposed 
by Gourgoulhon et al. p^ : 



W 



^ J CcntorNS 



(30) 

(31) 



W^ is the inertial velocity at the center of the star, and 
([30)) effectively separates the motion of the star relative 
to its center from its orbital motion. 

Note that Eq. ()29p is derived assuming the existence 
of an exact helicoidal Killing vector (for more details on 
the derivation of ([29]) from ([28]), read Teukolsky [Hi and 
Shibata [33|). This is, in general, not compatible with 
our choice of free variables in the XCTS equations. The 
error we introduce is most easily seen if we consider the 
evolution equation for the conformal factor. 



6 



(32) 



For Eq. ([M]) to be exact, we need 9tln0 — 0, while in 
the XCTS equations we assume that we are free to choose 
K = Q. As nothing guarantees that S/kP^ = — and 
in fact, we can check in practice that this term does not 
vanish — there is a contradiction within our equations^. 

Such approximations are inevitable, as there is no 
exact equilibrium solution to the binary problem. In 
practice, we will see that our numerical scheme is suf- 
ficiently accurate that they represent our main source 
of error. Better choices for A', or for our other free 
variables, might reduce these errors. However, within 
the quasiequilibrium formalism, we cannot hope to make 
them completely disappear. In fact, even though the 
contradiction here was shown using the hydrostatic con- 
ditions, a quasiequilibrium formulation creates very simi- 
lar problems in vacuum. (A discussion of deviations from 
quasiequilibrium in BBH binaries can be found in [27j . 
and the amplitude of the time derivative of the conformal 
factor observed there for irrotational binaries is similar 
to our results for BH-NS binaries.) 

Finally, to close our system of equations we need to 
choose an equation of state (EOS). Here, we will con- 
sider a poly tropic fluid, with polytropic index F = 2. 
The pressure P, energy and baryon density p and poi in- 
ternal energy epo, ^^nd enthalpy h then obey the following 



The most natural way to get rid of that contradiction would be 
to use equation I I32I I as the definition of K. The quantity dt \n<j> 
would then be a free variable, and could be set to 0. However, 
Pfeiffer showed [26l | that such a choice makes the operator of the 
XCTS system noninvertible. Alternatively, inserting 1321 in an 
iterative scheme driving dt In </) to seems to be unstable both 
for BBH [H and BH-NS binaries. 



relations: 
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h 
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= 1 + eH , 

Po 




P=(l + 


e)Po, 




epo = T 


P 



(33) 

(34) 
(35) 
(36) 



The method used, however, is independent of the EOS 
chosen — as long as, given /i, we can retrieve P, p, and 
po- Indeed, we only use the EOS to reconstruct the mat- 
ter quantities E, J, S, and po needed in Eqs. (fT9l) . (|20|l . 
(PT|) and (gni) from the enthalpy h. We use a T = 2 poly- 
trope as a reasonable first approximation to the nuclear 
equation of state, which will allow direct comparison with 
previous numerical results in Sec. lIVCI 



C. Boundary Conditions 

Building initial data for BH-NS binaries requires us to 
solve a set of elliptic equations: the constraints p^ . (|20p . 
and (|2ip and, in the case of irrotational binaries, an ad- 
ditional equation for the potential 'J, (P5|) . We thus have 
to provide boundary conditions at infinity and on the BH 
horizon for the XCTS variables 0, a(j), and /?', and on the 
surface of the NS for the potential '^ . 

At infinity (or, in practice, aX R = lO^^M, the outer 
boundary of our computational domain), we require a 
flat Minkowski metric in the inertial frame: 



/3o 

ad) 



0, 
1, 

: 1. 



(37) 
(38) 
(39) 



We excise the BH interior. Assuming that the BH is in 
equilibrium and that the excision surface is an apparent 
horizon leads to the set of conditions derived by Cook 
andPfeifi'er [l^: 



tViA^ch = -- ih'^V 
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'J 



Pa. 
P\ 
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if Si = a, 

(3' -Pa^s' = n^^xle'^^, 



(40) 

(41) 
(42) 



where s' 



^s* is the outward unit normal to the sur- 



face, /i*-' its 2-metric, x1 = Xi — ci are the Cartesian 
coordinates relative to its center, J is a projection of the 
extrinsic curvature on the excision surface defined in Eq. 
(28) of [23], and fi^^ is a free parameter determining the 
spin of the black hole. For a corotational BH, Vt^^ — 0, 
while the value required to obtain a nonspinning black 
hole is a priori unknown. A good first approximation, 
suggested in [231, is fi^^ = Vt, the orbital angular veloc- 
ity. This choice typically leaves the BH with a spin an 
order of magnitude lower than in a corotational binary. 
For better results, we follow the method introduced by 



Caudill et al. [3a| for BBH: we iterate over the value of 
QP^ to drive the BH spin to zero. This iterative method 
can be used to generate a BH of arbitrary spin. 

The last boundary condition required on the appar- 
ent horizon is only a gauge choice. However, that choice 
impacts the amplitude of the deviations from quasicqui- 
librium [27[. For conformally flat initial data, we will 
impose 



ds (acj)) = 0, 



(43) 



a choice that already gave good results for BBH bina- 
ries. We will discuss in Sec. IIVDI how this condition is 
modifled when we choose a different conformal metric. 

Finally, on the surface of the star, the boundary con- 
dition for ^ can be directly inferred from ([29]) : when the 
density tends towards 0, we are left with the equation 



%jdi^d.jP(i = 



h^nP'c^'' 



di Po- 



rn 



As Vpo should be along the normal to the surface of the 
star. (j44p is a boundary condition on the normal deriva- 
tive of ^. 



D. Orbital Angular Velocity 

In the construction of BH-NS initial data, the orbital 
angular velocity Q is, in general, a free parameter. In- 
deed, together with the initial radial velocity, it deter- 
mines the eccentricity and orbital phase of the orbit. 
Here, we consider binaries a few orbits before merger, 
where the trajectory is expected to be quasicircular. As 
a first approximation, we can require force balance at the 
center of the NS, as proposed by Taniguchi et al. [171 : 



V In /i = 0. 



(45) 



Force balance guarantees that the binary is initially in a 
circular orbit. As it neglects the infall velocity, it leads to 
a slightly eccentric orbit, but still constitutes a good first 
guess. Using Eq. (j^ . (|i5l) can be written as a condition 
on the lapse a and the Lorentz factors 7 and 70: 



Vln/i = V(ln^ I =0, 
ajj 



or, using the definitions ([M]) and (|26p . 

Vln(a2 -7,j/3*/3^) = -2Vln7. 



(46) 



(47) 



Effectively, this is a condition on the orbital angular ve- 
locity f2, if we remember that the shift is decomposed 
according to ([8]). Defining b to be the unit- vector along 
the axis passing through the centers of both compact ob- 
jects, we determine the angular velocity from 



6'V, ln(a2 _ ^^^p'^pj) ^ _26'V, ln7. 



(48) 



In theory, the angular velocity appears on both sides of 
the equation, but we only write explicitly the left-hand 
side, keeping 7 constant. We then check that fl converges 
when (j48p is inserted in our iterative solving procedure, 
described in Sec. IIIICI 

As we only solve (|47| along the direction b, we still 
have to impose force balance along the transverse direc- 
tions. To do so, we include a correction term when com- 
puting the enthalpy: if ho is the enthalpy computed from 
Eq. ([22|) , we use as the effective value of h 



h^ ho[l- (V^ln/io) • (r-CNs)], 



(49) 



where V^ = V — b(b • V) and cns is the location of the 
center of the NS. 

This choice drives the maximum of the enthalpy to- 
wards Cns- If the equilibrium was exact, V^ In/ig would 
vanish. For our quasiequilibrium binaries, its norm is less 
than 10-*^. 

An alternative method of imposing quasiequilibrium is 
to use the Komar mass Mk- If we have a timelike Killing 
vector, then Mk and A^adm, the ADM energy, should be 
equal. This condition is less convenient to impose during 
the solution, as global quantities like Mk and Madm 
cannot be reliably computed when the constraints are 
violated. However, we can use this equality as a test of 
our initial data, and verify that {Mk — Mabm) gets small 
as we converge. 

When we start applying the procedure described by 
Pfeiffer et al. [30] to reduce the eccentricity of the system, 
the situation is slightly different. We then prescribe the 
value of the orbital angular velocity as well as the initial 
radial velocity. Eq. P5|) is no longer useful. Instead, we 
adapt Eq. (|^^ so that it fixes the position of the star in 
all three spatial directions, replacing V_l by V. 

Note that if 9( is not an exact Killing vector, the equal- 
ity between Komar and ADM mass is lost. We can then 
use (Mif — Madm) only as an indicator of deviations from 
an exact equilibrium state. For low-eccentricity binaries 
with a nonzero infall velocity, those deviations are signif- 
icantly larger than when the angular velocity is fixed by 
Eq. (|48|) . and the infall velocity set to zero. 



E. Observing physical quantities 

We have just seen that, for quasiequilibrium configura- 
tions, computing the Komar mass and the ADM energy 
could be useful in finding the optimal angular velocity, or 
to ascertain how far from equilibrium our initial data are. 
To ensure that our initial configuration has the desired 
physical properties, a few additional quantities have to 
be computed. 

First, we want to be able to fix the mass of the compact 
objects. For a spinning BH, we define the irreducible 
mass Mglj, ADM energy in isolation M^^^, and spin 



parameter obh, 
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(50) 
(51) 
(52) 



where Jbh is the angular momentum of the BH. For the 
NS, we compute the baryon mass 
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NS 
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iiU'Ui 



rdV. 



(53) 



Here, 7 is the determinant of the conformal 3-metric 'jij. 
To check quasiequilibrium, we would like to know the 
ADM energy and the Komar mass of the system. Mea- 
suring the total angular momentum is also useful, mainly 
for comparisons with post-newtonian (FN) predictions or 
other numerical initial data. Those quantities are typi- 
cally defined as integrals on Soo, the sphere at infinity, 
which is not convenient for computations. Integrating by 
parts, we can transform these expressions into integrals 
on any sphere S enclosing all matter and singularities 
and, when needed, a volume integral on V, the region of 
our initial slice lying outside of S. Assuming conformal 
flatness, K = 0, and no constraint violations, this gives: 
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The decomposition into surface and volume integrals is 
not unique, but we found these expressions convenient, as 
the contribution of the volume terms decreases at least as 
l/r away from the center of mass, reducing our sensitivity 
to small numerical errors at spatial infinity. 

To make sure that the axis of rotation of the binary 
passes through the origin of our numerical grid, we also 
require that the ADM linear momentum vanishes. It is 



computed in a very similar way: 
1 



ADM 
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(57) 



and our solver moves the position of the BH center so 
that Padm is driven to zero. 

Finally, when discussing boundary conditions, we have 
seen that for irrotational binaries the correct value of the 
parameter ^Ibh is unknown. We thus need to find the 
value that makes the BH spin vanish. To compute the 
spin, we use approximate Killing vectors on the apparent 
horizon, following a method [39| similar to the work of 
Cook and Whiting [40|. 



F. Conversion to Physical Units 

In this paper, and in our numerical code, the system of 
units is based on the arbitrary choice of a unit mass: the 
ADM energy of the BH in isolation. Combined with the 
convention G = c = 1 , this choice is enough to determine 
all units of interest for BH-NS binaries. For applications, 
it is necessary to express results in astrophysical units. 
In this section, we give the conversion formulas. 

We first define the ADM mass of the neutron star 
Mj^g*^ as the ADM mass of an isolated NS of baryonic 
mass Afwd . The total ADM mass of the binary at infinite 



'NS- 



separation is then 



Mo = A/, 
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NS 
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ADM 
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and the mass ratio is defined as 

AfADM 

7? = — — — 

'^ Af ADM • 

Isolated neutron stars of given polytropic index are com- 
pletely described by their ADM mass and their compact- 
ness 

Ti.rADM 



c = 



'NS 



Rq 



(60) 



where Rq is the areal radius. Furthermore, stars of equal 
compactness but different masses are related by a simple 
scaling law. This can be seen by defining the length scale 

(61) 
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The whole problem is then invariant |41j under the trans- 
formation 

* (63) 
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and dimensionless quantity 
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In numerical simulations, we can thus retrieve all possible 
configurations by keeping only C and TZ as free param- 
eters, and choosing M^§^ = 1. Systems with different 
masses but the same neutron star compactness will obey 
the previous scaling, with 
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and i?* ly the value of i?poiy when Afgjj 
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We also define 



C(C) = 



i\ /ADM ' 

^"ns 



(67) 



a quantity which, for a given compactness, can easily be 
obtained from the solution of the Tolman-Oppenhcimer- 
Volkoff (TOV) equation. Then, if the baryon mass of the 
star is expressed in solar masses, 



M^g^mNsMg, 



the BH ADM energy wiU be 



M^DM ^ ^^j^sMq. 
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It is now straightforward to retrieve the meaning of our 
units of distance and time. A code distance d corresponds 
to the physical distance 

I, = <,|HCeU<ir5!:2E^xl.48k„. ,70) 
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while a code time t is equal to 
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Note, however, that for D to represent an actual physical 
distance, d has to be the proper separation 



ds. 



(72) 



and not the coordinate distance on our numerical grid. 

In our tests, we choose TZ = \ and k = 51.76, which 
gives C = 0.149 and C = 1.075. The conversion is thus 



D = rfr!!!^)x 1.79km 
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III. NUMERICAL METHODS 



Turn now to the numerical methods used to solve the 
initial data problem, and to the way the solver enforces 
simultaneously the various constraints on the system de- 
rived in Sec. [ni In this paper, we focus on the case of 
irrotational binaries with no initial radial velocity, even 



though the solver has also been used for single stars, coro- 
tational binaries, and infalling binaries. The chosen con- 
figuration is the most challenging of the four cases: the 
method for the other cases can be derived by omitting 
the irrelevant steps from what we present here. 

The core of the problem is the two sets of elliptic equa- 
tions, the XCTS system dH]) , dlQl) , and ^, and the ir- 
rotational condition on the potential 'i' (|29p. To solve 
these equations, we use the multidomain spectral elliptic 
solver (spells) developed by the Cornell-Caltcch collab- 
oration, as described by Pfeiffcr et al. 12411 . Improvements 
to SPELLS since the publication of [2J|, mainly the in- 
troduction of cylindrical subdomains, have increased its 
efficiency by about a factor of 3. The performance of 
the solver on distorted subdomains — such as a subdo- 
main with a boundary chosen to follow the surface of the 
neutron star — has also been improved, allowing us to 
solve the initial data problem in the presence of matter 
without Gibbs oscillations at the surface of the star. 

SPELLS has already been used successfully to solve the 
XCTS system for BBH binaries [H, ^. Here, when 
solving for the XCTS variables, we consider the matter 
terms as fixed, while in (P5)) . only the potential ^ is vari- 
able. We will detail in Sec. lIIICl how to combine the two 
groups of equations, as well as the additional conditions 
of force balance (|45p , vanishing ADM linear momentum 
and BH spin, and known BH and NS masses. But first we 
discuss some aspects of the solution of the elliptic equa- 
tions themselves: the numerical grid, and specifics of the 
irrotational potential equation. 



A. Domain Decomposition 

1. Numerical Grid 

The flexibility of the multidomain method used by 
SPELLS allows us to use relatively complex subdomain 
decompositions, adapting the numerical grid to the ge- 
ometry of the problem at hand. It also makes it possible 
to solve directly the whole XCTS system as a single set 
of nonlinear equations, without further decomposition of 
the XCTS variables, and using a relatively low number 
of grid points. 

For binaries in spells, we build the numerical grid 
from 14 subdomains, as follows (see Fig. [1]): 

• Around the BH, we use two concentric spherical 
shells and require their innermost boundary to be 
an apparent horizon. 

• The neighborhood of the NS is covered by an outer 
spherical shell with inner boundary mapped to the 
surface of the neutron star. This outer spherical 
shell touches an inner spherical shell which covers 
the whole neutron star, except a small region at 
the center. To avoid having to deal with regularity 
conditions at the center of a full sphere, the central 
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FIG. 1: Subdomain decomposition close to the compact ob- 
jects, in the equatorial plane. The apparent horizon of the BH 
(right) is an inner boundary of the numerical domain, while 
the surface of the NS (dashed line) is the boundary between 
the two spherical shells on the left. 



region is covered by a cube overlapping the inner 
spherical shell. 

• Three rectangular parallelepipeds cover the region 
surrounding the axis passing through the centers of 
the compact objects: one between the BH and the 
NS, and one on each side of the binary. 

• Five cylindrical shells around the same axis cover 
the intermediate field region. Their innermost 
boundary is, for three of them, within the paral- 
lelepipeds, and for the other two, within the outer 
shell surrounding each compact object. 

• The far-field region is covered by a spherical shell, 
with a l/r coordinate mapping allowing us to place 
the outer boundary at spatial infinity (or, in prac- 
tice, at i?= 10i°Af). 

At the second highest resolution, which we use as a 
reference to estimate the accuracy of the solution, the 
cube at the center of the star has 11x11x11 collocation 
points, the spherical shells around the compact objects 
have 19x18x36 points, the parallelepipeds 13x20x20 
points, the cylinders 14x15x13 (15 in the angular di- 
rection) or 14x15x20 (the higher resolution for the sub- 
domains closer to the compact objects), and the outer 
sphere 12x10x20. For comparison, the numerical grid 
used in [21| is built out of spherical shells with resolu- 
tion 41x33x32 or 49x37x36 around the black hole, and 
25x17x16 around the neutron star. 

To make convergence tests, we will need a single mea- 
sure of the resolution used. For a domain decomposition 
using subdomains with different basis functions and num- 
ber of collocation points, this definition is certainly not 



unique. Wc will use 



B. Irrotational flow 
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(75) 



where iV, is the number of collocation points in subdo- 
main i. For our second highest resolution, N^'^ = 44.0, 
while for Ref. [H ^^^^ > 78.7. 



2,. Surface Fitting 

Discontinuities in variables within a subdomain spoil 
spectral convergence. The surface of the star is a discon- 
tinuity, so we make it the boundary between two subdo- 
mains. (Note, however, that it is possible to reach a good 
level of precision — of the order of the error coming from 
deviations from quasiequilibrium — simply by including 
the surface in the interior of a thin spherical shell.) 

The surface of the star is approximated by an expan- 
sion in spherical harmonics. 
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where the center of the star, as defined in Eq. (|^5|) . is 
the origin of the spherical coordinates. To determine the 
coefficients c;,„, we solve the equation h{Rij,6i,(j)j) = 1 
along each collocation direction {6i, (j)j) of the numerical 
grid. Then, we project onto spherical harmonics the func- 
tion R{0, 0) defined by its values i?y in each collocation 
direction. 

To avoid Gibbs oscillations, we force the surface to be 
at the boundary between two spherical shells, Sq and 
Si. This is done by a coordinate transformation R — > 
R' fixing the radius of the common boundary between 
So and 5*1 to be the given function i?bound(^, 0)- This 
function is expanded in spherical harmonics, and will be 
equal to -Rsurf when the solver converges, as explained in 
Sec. IIII CI If 5*0 is defined in the original coordinates by 
Rq < R < R*, and Si hy R* < R < Ri, the map is, in 
So, 



R'{e,i 



-'Abound (^ 



Rn 



R*^Ro 



{R - Ro) + Ro, (77) 



while in Si we have 



K — rti 

The exact value of R* is not important, as long as 
Rq < R* < Ri. However, having R* ~ i?bound is usually 
convenient, as it leads to R ^ R'. 

The validity of the Yim expansion is evaluated by ob- 
serving the convergence of the coefficients cim as the res- 
olution increases. Results for a test irrotational binary 
are discussed in Sec. IIVBI 



Once the domain decomposition has been chosen, the 
XCTS equations can be solved without further modifica- 
tion. The irrotational equation (|29p . however, has spe- 
cific problems that require further attention. 

First, the coefficient of the leading order term — the 
Laplacian — vanishes on the surface of the star. As the 
equation is preconditioned by the inverse of a finite dif- 
ference approximation of the flat Laplacian, convergence 
will become extremely poor close to the surface, where 
(P^]) is very different from Laplace's equation. We thus 
change the preconditioning operator from an approxi- 
mation of — V^u to an approximation of —poV^u + u. 
The leading order term will then be properly represented 
within the star, while, when the density decreases, the 
operator becomes the identity and no preconditioning is 
done. 

Another problem is related to the inconsistencies in the 
quasiequilibrium formulation, already discussed in Sec. 
IIIBI Indeed, we know that, for a perfect equilibrium, 
Eq. (P^ will admit an infinite number of solutions (the 
potential is only defined up to a constant term). But, if 
we have instead a quasiequilibrium situation, Eq. (|29[) 
is not an exact representation of the continuity equation 
anymore. And nothing guarantees that a solution even 
exists. We found in practice that when using Eq. ([29]) 
as written, the convergence of the solver stops before we 
reach an acceptable precision. 

Different solutions to this problem were tried, involv- 
ing small modifications of Eq. ([29)1 . Here small means 
"at most of the order of the deviations from quasiequi- 
librium." The results presented here were obtained by 
replacing K in (j29p by the value required to ensure that 
(9f ln0 = using Eq. ([5^ . Of course, this does not solve 
the inconsistency — K \s, still set to in the XCTS equa- 
tions — but it guarantees that Eq. ([^5]) has a solution, 
allows the system to converge, and does not introduce 
any new source of error. 

Another method, mathematically less satisfactory but 
leading to equivalent results, is to allow for a small cor- 
rection in (j29p . for example, by adding the mean value 
of the potential, ^, to the boundary condition (pi)) and 
requiring that ^ is driven to zero (or, in practice, the 
small value required to counter the error coming from 
our choice of K) as we converge. 



C. Building quasiequilibrium binaries 

As discussed previously, knowing how to solve each set 
of elliptic equation is only part of the problem. Here, 
we outline how the solver links all of the requirements 
together and ensures convergence towards a solution rep- 
resenting the desired physical situation. 

At a fixed resolution, we solve according to the follow- 
ing algorithm: 
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1. Solve the XCTS system ([19]), (|20l), and ([2^, with 
fixed conformal matter quantities E, S and J*. The 
new value of the XCTS variables is determined by 
the relaxation formula u„ = (1 — A)w„_i + Xu*, 
where A is an arbitrary parameter (we typically use 
0.3) and u* the value of u found by solving the 
XCTS equations. In fact, knowing that we will 
use a relaxation formula, we do not even solve the 
equations exactly at each iteration; an approximate 
solution is good enough, and saves a lot of computer 
time. 

2. Impose symmetry across the equatorial plane (this 
step is not required, but we know that this symme- 
try should be respected, and enforcing it strictly 
accelerates convergence) . 

3. Evaluate the position of the surface of the star, 
-^surf' ^^'^ compare it to the evaluation made dur- 
ing the previous iteration, R"~^i- If both agree 
within a certain precision — we use the condition 

ll^sVf - R:J\\2 < 0.1||i?,Vf - i^boundlb, where 
-Rbound is the function used in the mapping ((77|) 
— modify the numerical grid by setting i?bound = 



K 



surf ■ 



4. Compute the ADM linear momentum P 



and 
compare it to the value computed during the pre- 



ADM' 



vious iteration, 



pri-l 
^ADM- 



If IIP 



ADM 



PadmII < 



0.1 X IIPadmII' ™ove the center of the BH. The 
change in the position of the center, (5c, is chosen 
so that, if the system was Newtonian, the total lin- 
ear momentum would vanish: (5c x J7 = Padm- 
We also change the radius of the excision surface 
(the inner boundary of the shells around the BH) 
to drive M^^^ to its desired value. 

5. Solve Eq. (P5|) to find the new angular velocity. 

6. Get the spin of the BH, and change the parameter 
57^^ in the boundary condition P^ to drive the 
spin to — or any other desired value, if the BH 
is not irrotational. The new value of fi^^ is chosen 
by linear interpolation, using the last two values of 
the spin. 

7. Determine the constant in the Euler first integral 
(P^ so that the baryon mass of the NS ([55)) is set 
to its target value. 

8. Apply correction ([49]) to the value of the enthalpy. 



9. Solve the irrotational equation ()29p for "$. The new 
value of V? is determined using the same relaxation 
formula as for the XCTS variables. 

10. If the desired precision has not been reached, go 
back to 1. 

From this description, it is clear that the accuracy of 
the results depends on the convergence of the many pa- 
rameters updated during the iterative procedure. We will 



TABLE I: Domain decomposition for a single TOV star. For 
spherical shells, the three numbers denote the resolution in 
radial, polar and azimuthal directions. 





Central Cube 


Inner Shells 


Outer Shell 


RO 


7x7x7 


7x6x12 


8x6x12 


Rl 


9x9x9 


10x9x18 


9x7x14 


R2 


11x11x11 


13x12x24 


10x8x16 


R3 


13x13x13 


16x15x30 


11x9x18 



discuss in Sec. IIVBI various tests verifying that they all 
reach an acceptable precision. 



IV. TESTS AND RESULTS 

As mentioned earlier, the main motivation to build a 
code generating BH-NS initial data using a multidomain 
spectral method is the possibility of rapidly reaching high 
levels of precision. As an example, we will focus on a se- 
quence of irrotational, equal-mass BH-NS binaries. In 
Sec. IIVBI we show through convergence tests that, over 
a large range of initial separations likely to be chosen as 
starting points for future evolutions, we can construct 
initial data with enough precision to resolve deviations 
from quasiequilibrium. Trying to reach higher precision, 
even if mathematically possible, would be of little inter- 
est: the additional information would not be physically 
meaningful. 

We then turn, in Sec. llVCi to another interesting test 
of our results: comparing them to a similar sequence gen- 
erated by Taniguchi et al. [21|, and to predictions from 
the 3PN approximations computed by Blanchet [42|, as 
well as Mora and Will [43. With accurate estimates of 
our errors, we discuss how far deviations of the numerical 
results from the 3PN approximations can be trusted, and 
their potential interpretation. 

Finally, we end this section with a discussion of two 
different types of initial configurations: binaries built us- 
ing a modified Kcrr-Schild conformal metric to construct 
systems with a spinning BH, and binaries with an initial 
radial velocity, which can be used to generate systems 
with low-eccentricity orbits. 



A. TOV Star 

Before tackling binaries, we test our algorithm on an 
isolated, nonrotating NS. This effectively means that only 
steps 1, 2, 3, and 7 of our solution procedure are not 
trivial. Although the position of the surface is known 
analytically, for the purpose of this test we rely on the 
iterative surface fitting method to find it. An "exact" 
solution is easily computed by direct integration of the 
TOV equations. We compared the central density, ADM 
mass, Komar mass, and central lapse: all converge expo- 
nentially with resolution. 
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FIG. 2: Error in the energy density for an isolated NS as a 
function of the isotropic radius Riso- The reference configu- 
ration is obtained by numerical integration of the TOV equa- 
tions. The spikes in the error function are due to a change in 

the sign of (pnum - Pan) 



Figure [2] shows the difference between the exact and 
computed density profiles. We can see that the spectral 
convergence of the error holds at all radii. 

For this simple case, the domain decomposition con- 
sists of just a cube covering the center of the star, two 
spherical shells whose common boundary matches the 
surface, and a third shell with an l/r mapping extending 
toR= IQi^M. The resolutions RO to R3 used in the test 
are described in Table [J 



B. Irrotational binaries 

To test the performance of our solver for binary sys- 
tems, we use the iterative method from Sec. 1111 CI to 
construct a sequence of equal-mass, irrotational binaries. 
The NS has an ADM mass in isolation of 1 (in code units: 
see Sec. Ill FI for a conversion in astrophysical units), and 
a parameter k = 51.76, leading to a compaction simi- 
lar to that used in Ref. [2l|, Table IV. Our results are 
detailed in Table M 

We look at three different sources of error. 

• The iterative procedure. To estimate that error, 
we study the convergence, at fixed resolution, of all 
the parameters changing between iterations. 

• Truncation errors. We observe the convergence of 
the solution with the number of collocation points 
by solving each configuration at four different reso- 
lutions, RO to R3, as detailed in Table Hill The sec- 
ond highest resolution, R2, is our standard numer- 
ical grid, as defined in Sec. 1111 A H and the highest 
resolution, R3, is used as an approximation of the 
exact solution. 



TABLE II: Sequence of irrotational, equal-mass BH-NS bina- 
ries. We give here the coordinate distance between the centers 
of the two compact objects d, the orbital angular velocity fl, 
the binding energy Ei,, the angular momentum J, and the 
difference between Komar and ADM energies. 
For the four closest configurations, marked by an asterisk, 
the numerical error estimated from the convergence of energy 
measurements is larger than the deviations from quasiequilib- 
rium, approximated by 5M = Mk — Eadm, so that 5M might 
not be resolved. The error in Eb reaches about 5 x 10~^ 
at d/Mo = 8.406, an order of magnitude larger than at 
d/Mo = 9.007. 



d 
Mo 


nMo 


Mo 


J 


\Eadm-Mk\ 
Mo 


18.505 


0.01169 


-6.1490 X 


10"=^ 


1.19460 


7.7 X 10-^ 


16.506 


0.01377 


-6.8103 X 


10"=^ 


1.14241 


9.5 X 10-^ 


14.506 


0.01653 


-7.6289 X 


10-3 


1.08815 


1.4 X 10-® 


12.506 


0.02037 


-8.6634 X 


10"=^ 


1.03177 


2.0 X 10-® 


11.506 


0.02288 


-9.2879 X 


10-=* 


1.00284 


2.6 X 10-® 


10.507 


0.02596 


-1.0002 X 


10-2 


0.97353 


3.4 X 10-® 


9.507 


0.02981 


-1.0821 X 


10-2 


0.94408 


4.4 X 10-® 


9.257 


0.03092 


-1.1043 X 


10-2 


0.93675 


5.0 X 10-® 


9.007 


0.03211 


-1.1273 X 


10-2 


0.92947 


5.4 X 10-® 


8.857 


0.03285 


-1.1416 X 


10-2 


0.92514 


5.6 X 10-®* 


8.757 


0.03337 


-1.1509 X 


10-2 


0.92225 


6.0 X 10-®* 


8.557 


0.03445 


-1.1706 X 


10-2 


0.91656 


6.6 X 10-®* 


8.406 


0.03530 


-1.1853 X 


10-2 


0.91237 


8.2 X 10-®* 



• Deviations from equilibrium. We know that the 
quasiequilibrium formalism contains intrinsic con- 
tradictions. A useful estimate of the error thus 
created is the difference between the Komar and 
ADM energies. In the presence of an exact time- 
like Killing vector, both would be equal, but here 
the difference can be seen as an indication of how 
far from equilibrium we are. 

All the graphs presented in this section correspond to 
a binary with rescaled coordinate separation cI/Mq = 
11.507. A summary of our results for the whole sequence 
is in Table [III Typically, the numerical error rises as the 
separation decreases. The difference between Komar and 
ADM mass can be resolved up to d/Mo = 9. Numerical 
errors then start to increase rapidly to reach, for our clos- 
est binary, values around 5 x 10"^. By that point, the 
solver does not converge at resolution R3 anymore, and 
we thus use Rl as our reference and R2 as an estimate 
of the exact solution. 



Convergence of the iterative 



To verify the convergence at fixed resolution, observe 
Figs. [3] and [H The iterative procedure converges if all 
parameters modified within one step converge, while the 
residuals from the two elliptic solves (i.e., the constraint 
violations and the deviations of the fluid from an irro- 
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TABLE III: Domain decomposition for binary systems. A 
description of the different subdomains can be found in Sec. 
IIII A II Tfie tfiree numbers denote tfie resolution in radial, 
polar, and azimuthal directions for spherical shells, and in ra- 
dial, polar, and axial directions for the cylinders. The cylin- 
ders have two different resolutions (HR/LR), the highest be- 
ing used for the two subdomains directly surrounding one of 
the compact object. Finally, for the parallelepipeds, the first 
number corresponds to the resolution along the axis passing 
through the centers of both compact objects. 





Cube 


Inn. Shells 


Out. SheU 


Par all. 


Cyl. (HR/LR). 


RO 


9x9x9 


13x12x24 


8x6x12 


9x12x12 


10x9x12/9 


Rl 


10x10x10 


16x15x30 


10x8x16 


11x16x16 


12x12x16/11 


R2 


11x11x11 


19x18x36 


12x10x20 


13x20x20 


14x15x20/13 


R3 


12x12x12 


22x21x42 


14x12x24 


15x24x24 


16x18x24/15 




50 100 

Step 



150 



tational configuration) vanish. In Fig. [3l we show the 
evolution of three of these parameters while iterating at 
our lowest resolution RO: the angular velocity fl, derived 
from the Eq. (|48|) , the constant in the Eulcr first integral 
([22|) . which controls the mass of the NS, and the areal 
mass of the BH, controlled by the radius of the excision 
surface. The difference between the parameter at a given 
step and its final value at the highest resolution is shown. 
We see that, even though the resolution is low, all param- 
eters converge to a relative precision below 10~^. At the 
reference resolution R2, the relative precision is better 
than IQ-'^. 

In addition to the overall convergence. Fig. [3] shows 
abrupt changes, especially in the evolution of the BH 
mass. These can easily be understood if we remember 
how the mass of the BH is fixed: the radius of the ex- 
cision boundary is modified whenever the linear ADM 
momentum converges. We then change our numerical 
grid and the location of the apparent horizon. Every 
time we regrid, the BH mass will at first be very close to 
its desired value, then reach a new equilibrium when the 
system adapts to its new boundary condition. The mass 
just before regridding — when the error is maximal — is 
thus the best estimate of our precision. 

We also monitor the evolution of a number of quan- 
tities that should tend towards zero as the system con- 
verges: the total linear momentum (to ensure that the 
axis of rotation passes through the origin of our coordi- 
nate system), the BH spin (as we want irrotational bina- 
ries), the quantity V_l In ft in Eq. (|1^ . and the L2 norm 
of modes violating the equatorial symmetry (before we 
manually impose it). The last converges quickly to rela- 
tive precisions of order 10~^, and down to about 10~^° 
at resolution R2, while the behavior of Padm and Jbh is 
shown on Figured As in Fig. [31 we plot the evolution at 
our lowest resolution, RO. We observe rapid convergence, 
with once more some oscillations due to the occasional 
modification of the numerical grid. At the reference res- 
olution R2, both Padm and Jbh vanish to a precision 
better than 10^^. From Figs. [3] and 21 we can thus safely 



FIG. 3: Convergence of the angular velocity, the Euler con- 
stant (which controls the mass of the star) and the mass of the 
BH while iterating at the lowest resolution RO for an equal- 
mass binary with initial separation d/Mo = 11.507. The val- 
ues plotted are the differences from the final results at the 
highest resolution R3. One step is defined as a passage from 
point 1 to point 10 in the iterative procedure described in Sec. 

lilTcl 




150 



Step 



FIG. 4: Convergence of the spin of the BH Jbh and the total 
linear momentum Padm at our lowest resolution RO for an 
equal-mass binary with initial separation d/Mo — 11.507. 



consider that the iterative method detailed in Sec. IIII CI 
does indeed converge at fixed resolution. 

The last parameter, V±\nh (not plotted), does not 
however completely vanish, even at our highest resolu- 
tion. In fact, it converges rapidly towards a fixed, small 
value of order 10~^. This is most likely because the equi- 
librium is not perfect — and, indeed, when the deviations 
from exact equilibrium increase, so does the final value 
of Vi In/i. 
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FIG. 5: Convergence of the Hamiltonian and momentum con- 
straints with resolution for an equal-mass binary at a separa- 
tion d/Mo = 11.507. 
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FIG. 6: Convergence of the surface fitting method measured 
as the evolution of the error in the coefficients of the expan- 
sion of -Rsurf in spherical harmonics, computed here as the 
difference with our results at our highest resolution. 



2. Spectral convergence of the solution 

Having established that the iterative procedure works 
as intended, we turn to an estimate of the precision of the 
initial data obtained, that is, the differences between the 
solutions at different resolutions. As we use a spectral 
representation, we expect exponential convergence of all 
variables. We report the convergence of the constraint 
violations, the performance of the surface fitting method, 
and the convergence of a set of measured global quantities 
(A/adm, t^ADM, MK^ and the position of the BH center 
cbh)- 

Fig.[5]shows the residual of the elliptic equations corre- 
sponding to the Hamiltonian and momentum constraints. 
At the end of an elliptic solve at any given resolution, it 
should vanish at all collocation points. In order to obtain 
a meaningful estimate of the error, we thus evaluate the 
residual on the numerical grid corresponding to the next 
higher resolution. The exponential convergence is clearly 
seen, and we can deduce from Fig.[5]that the norm of the 
constraints at resolution R2 is around 10~^. 

The performance of the surface fitting method can be 
evaluated from Fig. [51 where we show the convergence 
of the surface at different resolutions. The error is es- 
timated by the L2 norm of the difference between the 
coefficients of the expansion in spherical harmonics (j76p 
at the current resolution and their final values at our 
highest resolution. The exponential convergence allows 
us to easily estimate the error in the position of the sur- 
face. For this configuration the position of the surface is 
known within better than 10~^ code units. For highly 
distorted stars however, this error becomes significant, 
and provides the easiest way to check during the compu- 
tation whether the angular resolution is high enough or 
not. 
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FIG. 7: Convergence of the ADM energy, ADM angular 
momentum, difference between Komar and ADM energies, 
and position of the center of the BH with resolution for an 
equal-mass binary at a separation d/Mo — 11.507. We plot 
AEb = {Et, - E;)/Mo, AJ = (J - J*)/Mo^ A(£; - K) = 
{SM — SM*)/Mo, and Ac = ||cbh — Cbh)I|2, where the refer- 
ence results E^, J* , SM* and Cgfj are those at resolution R3 
(N^/^ = 51.5). 

The difference SM between AIk and -Eadm, an indication of 
how close to equilibrium the system is, reaches 2.6 x 10^^ at 
the highest resolution. This is significantly larger than the 
estimated error in either -Eadm, or SM shown in the figure. 



Finally, in Fig. [7] we show the convergence of the mea- 
sured ADM energy and angular momentum with resolu- 
tion. For both quantities, the reference for comparison is 
the value measured at the highest resolution R3. We see 
good convergence over 2 orders of magnitude. Similar 
figures can be obtained for different binary separations 
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— though as discussed earlier, our abihty to solve accu- 
rately at high resolution decreases when the star becomes 
too distorted. 



3. Deviations from equilibrium, 

Also in Fig. [71 we plot the convergence of the position 
of the BH center, which confirms that the center of the 
numerical grid is indeed the center of rotation of the sys- 
tem, and the convergence of the difference between the 
ADM and Komar energies 5M , a measure of the devia- 
tion from quasiequilibrium. Wc sec that this difference is 
resolved to a very high precision, much lower than its ac- 
tual value of 2.6 x 10^®. As the ADM energy itself is also 
resolved to a precision significantly better than 10~^, we 
see that our main sources of imprecision are the inconsis- 
tencies inherent in the quasiequilibrium approximation. 

Both the numerical errors and the deviations from 
quasiequilibrium increase as the separation decreases, 
but, as the star approaches its mass-shedding limit, the 
numerical error increases much more rapidly. As pre- 
viously mentioned, they are roughly comparable for a 
rcscalcd coordinate separation cL/Mq = 9. The decrease 
in performance at lower separations is not, however, a se- 
rious problem. By that point, the radial velocity of any 
real binary will already be significant, and so would other 
deviations from the idealized quasiequilibrium state. Any 
evolution looking for such levels of precision should prob- 
ably start from a larger separation. Results at small co- 
ordinate separations are, however, interesting for more 
qualitative predictions. For example. Taniguchi ct al. 
[2l| use them to determine which configurations are likely 
to reach the innermost stable orbit before the star gets 
disrupted. We will thus keep them as useful approxima- 
tions, without expecting the same precision as for more 
widely separated objects. 



C. Comparison with previous results 

As a last test of our code, we compare the initial data 
generated using the iterative method described in Sec. 
unci to 3PN approximations and previous numerical re- 
sults. For these comparisons, we use the sequence of 
equal-mass, irrotational binaries detailed in Table HIl The 
3PN values were obtained in the point-mass, circular or- 
bit approximation by Blanchet [4^. We also use results 
from Mora and Will [4J| to take into account eccentric- 
ity and finite size effects. For the numerical comparison, 
we use the data from Table 1 V of Taniguchi ct al. [2l| . 
These last results are given to 3 significant digits, the 
actual precision being unknown to us. Their error in the 
quasiequilibrium condition — our sole basis for compar- 
ison — is, at most separations, around an order of mag- 
nitude higher than what we observe in our initial data. 
This error is, however, small enough to allow comparisons 
of both numerical results with the 3PN approximations. 



TABLE IV: Choice of 3PN models used as references. The 



iccentriciti 


^ e is defined as in 


a, Eq. (2 


3) 






Source 


Finite size 


e 


Orbital pos. 


3PN-B 


Blanchet [42j 


No 





— 


3PN-M0 


Mora and Will [43] 


Yes 





— 


3PN-MP 


Mora and Will [4^ 


Yes 


0.01 


Pericenter 


3PN-MA 


Mora and Will [43] 


Yes 


0.01 


Apocenter 



Four different models are compared. The first cor- 
responds to the results of Blanchet [43, where the 
orbits are circular and the compact objects are modeled 
as point masses. The second adds finite size effects 
to the model. Most corrections made by Mora and 
Will [43[ to the point-mass model vanish in the case of 
an irrotational binary, and only the tidal effects add a 
significant contribution. We compute them according 
to Eq. (3.6a) of their work. The last two models 
include some eccentricity. The exact eccentricity of our 
initial data is, in general, unknown. However, we can 
get reasonable estimates from evolutions starting at 
separation d/A/o = 12. We will give the 3PN results 
for binaries with an eccentricity e = 0.01. At a given 
eccentricity, the binding energy and ADM momentum 
reach extrcma at the pericenter and the apocenter. We 
thus present the 3PN results at those two points, giving 
an order of magnitude estimate of the influence of the 
eccentricity. A summary of the parameters chosen for 
the four models is given in Table [TVl 

In Fig. [51 we show results for the binding energy for 
various binary separations, where both numerical simu- 
lations seem to be in good agreement. For our results, 
the precision reached is good enough to measure devi- 
ations from the 3PN predictions neglecting eccentricity. 
We observe differences of order 10"^ for configurations 
where our expected precision is about an order of magni- 
tude better. In Fig. [51 we show the deviations from the 
simplest 3PN model (3PN-B) over a large range of sepa- 
rations. The behavior at small separation is not resolved 
well enough to note anything other than the divergence of 
the numerical and 3PN predictions when the star reaches 
its disruption point. But for most of the sequence, we 
observe that the numerical results are clearly below the 
3PN predictions, the difference between the two results 
decreasing at the largest separations. 

It is also easy to see that tidal effects cannot explain 
these results. They contribute at the same order of mag- 
nitude, but tend to increase the energy of the system. 
However, Fig. [S] shows that our results are still compati- 
ble with the 3PN predictions if we include the influence 
of eccentricity. Indeed, its effects can decrease the energy 
of the system if we are closer to the apocenter than the 
pericenter — and, in fact, we know from short evolutions 
that this is the case for our initial data (see Table IVlI)) . 

A similar comparison can be made using the total an- 
gular momentum Jadm, as shown in Fig. 1101 The agree- 
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FIG. 8; Binding energies of equal-mass binaries for initial 
data from our solver, from Taniguchi et al. [2l[, and from 
3PN predictions for model 3PN-M0 (see Table ITV)). 
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FIG. 10: Angular momentum of equal-mass binaries for initial 
data generated by our solver, initial data from Taniguchi et 
al. [2]| . and 3PN predictions. We see that both numerical 
results are in very good agreement, even when they begin to 
diverge from the 3PN models of circular orbits. 



3PN predictions if we assume a small eccentricity and an 
initial state closer to the apoccnter than to the pericen- 
ter. 

Overall, these results show that our precision is good 
enough to resolve deviations from the point-mass, circu- 
lar orbit 3PN predictions. We have thus the potential to 
study the main effects contributing to these deviations in 
irrotational BH-NS binaries: tidal effects in tlie neutron 
star, influence of the eccentricity of the orbit, and spuri- 
ous gravitational effects due to the inconsistency of the 
quasiequilibrium formulation. 



FIG. 9: Difference between our results for the binding energy 
of a sequence of quasiequilibrium equal-masses binaries, and 
the 3PN predictions from model 3PN-B. The errors repre- 
sented here come from the difference between the ADM and 
the Komar energies, except for the 3 closest binaries, for which 
the numerical error can no longer be neglected. 
We also represent the influence of tidal effects (from model 
3PN-M0) and eccentricity (model 3PN-MA). Any binary with 
an eccentricity e — 0.01, initially closer to its apocenter than 
to its pericenter, should have a binding energy between the 
results from models 3PN-M0 and 3PN-MA. Model 3PN-MP, 
representing an eccentric binary at its pericenter, is not plot- 
ted here, but predicts even higher energies than model 3PN- 
MO. 



ment between both numerical calculations is clearly visi- 
ble, even in the regime where they deviate from the 3PN 
models of circular orbits. This should not be surprising, 
as both sets of numerical results use essentially the same 
formulation of the problem. As was the case for the en- 
ergy, results for Jadm can only be reconciled with the 



D. Spinning black holes 

For nonspinning or slowly spinning BHs, the confor- 
mally flat metric we have used until now performs ex- 
tremely well. However, if we want to generate rotating 
BHs, being able to use a different conformal metric is 
critical. The natural choice for such BHs would be a 
Kerr-Schild conformal metric. Unfortunately, early re- 
sults from Taniguchi et al. [l^ have shown that this 
leads to strong deviations from quasiequilibrium: the dif- 
ference between the Komar and ADM energies in [l^ is 
of tlie order of the binding energy. Here, we describe the 
use of a modiflcation of the Kerr-Schild metric, already 
applied to BBH by Lovelace et al. [33j . 

We deflne 7,^'^(aBH, vbh) and if^'^(aBH, vbh) as the 
3-metric and trace of the extrinsic curvature of a black 
hole with spin parameter cbh and boost velocity vbh, 
written in Kerr-Schild coordinates. Then, we choose the 
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free parameters of the XCTS equations as follows: 

l^o - '5,, + [7,f^(aBH,VBH)-'5.,]e-('-^/'"'', (79) 



vbh = rz X Cbh 



(80) 
(81) 



where ri is the coordinate distance to the center of the 
BH Cbh , and the width w is chosen as half the coordinate 
distance between the two compact objects. This choice 
ensures that close to the BH, the metric is nearly J^^ , 
while away from the hole, we recover conformal flatness 
and maximal slicing. The introduction of the exponential 
damping e"^''^/"'^ is the most important difference be- 
tween the choices of conformal metric and extrinsic cur- 
vature in [33[ and [l9| . That change is indeed necessary 
to avoid large deviations from equilibrium. 

To take advantage of the similarities between this ini- 
tial conflguration and a Kerr black hole, we also change 
the boundary condition imposed on the lapse. If the 
lapse of an isolated Kerr-Schild BH is Q;''^'^(aBH, vbh), 
our boundary condition on the excision surface will be 



,KS 



(aBH,VBH)e 



-(n/w) 



(82) 



instead of (gS]). To get as close as we can to a Kcrr- 
Schild BH, we modify the shape of the excision surface. 
The subdomain containing the apparent horizon is now a 
spherical shell in coordinates (tk, 0,4')- The Kerr radius 
tk is defined as the largest positive root of the equation 
rj^ — rj^{r^ — a^) — {a. ■ r)^ =0, where r is the coordinate 
distance to the center of the BH, and a = J /AI^^^ is the 
spin parameter. The excision surface is then the oblate 
surface rn = constant, and we choose the constant so 
that A/ADM = 1 

Once these choices have been made, no further modi- 
fications of our numerical methods are required. We test 
the performance of these new data sets on two types of 
BH-NS binaries. First, we consider configurations with 
a nonspinning BH, which allows direct comparison with 
the conformally flat initial data. Then, we move to BHs 
with a spin Jbh = 0.5(Afg^'^)^, with the direction of 
Jbh opposite to the orbital angular momentum, and ver- 
ify that comparable results can be obtained. Tables IVl 
and IVII summarize the properties of the resulting bina- 
ries. Different spins aligned with the rotation axis can be 
obtained using the same method. We tested our proce- 
dure up to spins of 0.9, and note that, as for BBH initial 
data [23], the deviations from quasiequilibrium tend to 
increase with the spin of the BH (the difference between 
Komar and ADM mass reaches about 10% of the binding 
energy for a spin of 0.9). The choices of the conformal 
metric and the lapse boundary condition seem to have 
a major influence on the amplitude of these deviations. 
Better choices will probably help reduce the deviations 
observed for rapidly rotating BHs. 

We first note that, for equivalent resolutions, the new 
configurations are less precise by typically an order of 
magnitude. Also, as we want the width w to be large 



TABLE V: Same as Table HH but for BH-NS binaries built 
with a modified Kerr-Schild conformal metric, as described in 
Sec. ITVDl The spin of the BH is still 0. 



d 

Mo 



18.489 
16.990 
15.490 
13.991 
12.491 
11.492 
10.493 



fiM, 







0.01171 
0.01321 
0.01507 
0.01741 
0.02042 
0.02295 
0.02605 



A/f, 



-6.15 X 10"=* 
-6.64 X 10-3 
-7.20 X 10-3 
-7.87 X 10-3 
-8.67 X 10-3 
-9.30 X 10-3 
-1.00 X 10'^ 



ALj 



1.195 
1.155 
1.116 
1.074 
1.032 
1.003 
0.974 



Mo 



8.4 X lO"'' 
8.8 X 10-*^ 
9.3 X 10"^ 

1.0 X 10-^ 

1.1 X lO-'"^ 
1.1 X 10"^ 
1.3 X 10"^ 



TABLE VL Same as Table |Vl but the BH now has a spin 
Jbh = -0.5. 



d 
Mo 


QMo 


Ma 


M- 


Badm-A/a' 
Ma 


18.368 


0.01182 


-6.03 X 10^3 


1.081 


2.9 X 10"^ 


16.881 


0.01335 


-6.50 X 10^3 


1.043 


3.4 X 10^^ 


15.395 


0.01523 


-7.04 X 10-3 


1.004 


4.0 X lO"'"^ 


13.908 


0.01759 


-7.68 X 10-3 


0.964 


4.7 X 10"^ 


12.422 


0.02065 


-8.43 X 10-3 


0.924 


5.6 X 10-^ 


11.431 


0.02321 


-9.01 X 10-3 


0.896 


6.2 X lO-"^ 


10.441 


0.02636 


-9.67 X 10^3 


0.869 


6.8 X 10^^ 



compared to the radius of the apparent horizon, we 
should avoid close binaries. Deviations from quasiequi- 
librium are also significantly larger, but not nearly as 
much as in [l9[ , where an unmodified Kerr-Schild back- 
ground was used. In [l^, the difference between Komar 
and ADM energies was of the order of the binding en- 
ergy, while here it is only about 0.15% of that value. Di- 
rect comparison between our results for a flat conformal 
metric and Table |V] also shows that both sets of initial 
configurations are in agreement. 

As long as the BH is not rotating, the new conformal 
metric does not lead to any noticeable advantage over the 
conformally flat background — though the initial burst of 
gravitational radiation might end up being smaller. For 
rotating BHs, however, a conformally fiat metric is no 
longer appropriate, while a modified Kerr-Schild metric 
allows us to solve the initial data problem. Deviations 
from quasiequilibrium will increase once more, but for a 
BH spin Jbh = —0.5, we can still solve for the binding 
energy within a fraction of a percent. The norm of the 
constraints is also below 5 x 10-^ for our closest binaries, 
making these initial configurations perfectly suitable for 
future evolutions. 

Our ability to reach high accuracy at a relatively low 
resolution is particularly important for the construction 
of these spinning configurations. Indeed, the slower con- 
vergence rate makes it significantly harder to obtain use- 
ful initial data. Moreover, as the geometry around each 
compact object become less and less spherical, being able 
to easily adapt our numerical grid becomes even more 
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TABLE VII: Orbital parameters of three irrotational BH-NS 
binaries, after 0, 1, and 2 steps of the iterative procedure 
designed to reduce the eccentricity of their orbits. The initial 
radial velocity of an observer comoving with the NS is Vr ~ 
aorfo/2, the eccentricity is measured from the parameters of 
the fit (|83|) according to e = B/codo, and the orbital phase (p 
is at pericenter and n at apocenter. 





Vr 


QMo 


e 


<t>/^ 


\Eadm-Mk\ 
M„ 


StepO 
Stepl 
Step2 



-9.36(-4) 
-7.20(-4) 


0.02157 
0.02161 
0.02165 


1.0 X 10"^ 

4.4 X 10"^ 

6.5 X 10~* 


0.68 
1.18 
1.59 


2.3 X 10"'^ 

2.8 X 10"* 

2.9 X 10"* 
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FIG. 11: Evolution of equal-mass binaries after 0, 1 and 2 
steps of our iterative method reducing the eccentricity. We 
plot the time derivative of the coordinate separation between 
the BH and the NS, 2iJr = d. 



necessary. These first results show that the construction 
of spinning BHs is perfectly possible without much mod- 
ification of our basic formalism — and improvements in 
the choice of the conformal metric and/or the excision 
boundary conditions might further improve the quality 
of such initial configurations. 



E. Low-eccentricity binaries 

The initial configurations discussed in this paper cor- 
respond to binaries only a few orbits away from merger. 
Such systems are expected to have nearly circular orbits 
as, because of gravitational wave emission, the eccentric- 
ity decreases as a power law of the distance between the 
44| . The influence of the eccentricity on observ- 



objects 

able quantities such as the gravitational waveform can be 

significant. For instance, it is one of the dominant effects 



limiting the comparison between high-accuracy BBH evo- 
lutions and post-Newtonian expansions presented in [4^ , 
even though the initial eccentricity of their binary is lower 
than 6 x 10"^ 

Evolutions of BH-NS systems are far from being as 
precise. But the force balance condition we used until 
now leaves the binaries with eccentricities of order 0.01 
— enough to be noticeable in evolutions. We thus want 
to decrease the eccentricity of the initial data so that 
its influence on the orbit is at most of the order of the 
precision of the evolution code. 

Here, we show that the iterative method already used 
to reduce the eccentricity of BBH [30] can be applied suc- 
cessfully to BH-NS binaries. For all evolutions described 
in this section, we used the mixed finite difference- 
spectral code described in Duez et al. |4q . 

The eccentricity and orbital phase of our binaries are 
determined by the choice of orbital angular velocity f2 
and infall velocity dor. Until now, we have been deter- 
mining n through Eq. (|l5)) . choosing do = 0. Now, we 
will use such configurations as a first approximation to 
the low-eccentricity solution, and try to determine from 
its evolution better values of CI and do- 

To do so, we record the coordinate separation between 
the center of the compact objects, d, and fit its time 
derivative by the formula 



d = Aq + Alt + B sin {ujt 



(83) 



where the parameters ^Oi Ai, B, w, and arc all deter- 
mined by the fit. For a Keplerian orbit, we would have 
Aq = Al = 0, and an eccentricity e = B/cudo, where 
do = d{t = 0). We use this definition of e as an approxi- 
mation of the eccentricity of the system. As in [30| , we 
then choose the corrections to fl and do so that a Ke- 
plerian orbit with the same parameters d, uj, 0, and B 
would become circular: 



Sdo = 

sn = 



Bsincj) 

do 
Bu! cos 4 

zdo* *o 



(84) 
(85) 



For the fit (|85|) to be accurate, we need to evolve the 
binaries for at least one and a half orbits. Furthermore, 
as the initial spurious burst of gravitational radiation in 
the data disturbs the early motion of the binary, we also 
exclude points at i < lOOM from the fit. 

As a first example, we consider a binary at initial co- 
ordinate separation d/Mo = 12.0, and evolve it using the 
fully relativistic numerical code described in [46|. From 
this evolution, we determine that the eccentricity of the 
initial data constructed by requiring force balance ((45|) 
and dor = is of order e = 0.01. We then go twice 
through the iterative method we just described. The 
orbital parameters of the three binaries we evolved are 
listed in Table IVIII while in Fig. (TTJ we show the time 
derivative of the coordinate separation, d. Two itera- 
tions reduce the eccentricity by about an order of magni- 
tude. Decreasing the eccentricity further would demand 
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evolutions at a higher resolution, increasing the computa- 
tional cost, but does not in principle involve any new dif- 
ficulties. We also find that the difference between ADM 
energy and Koniar mass increases by about 2 orders of 
magnitude during eccentricity removal (see Table IVIIj) . 



V. DISCUSSION 

In this paper, we presented a new method for the con- 
struction of initial data for BH-NS binaries, based on the 
multidomain spectral elliptic solver spells [2J]. The 
flexibility of the multidomain spectral methods allows the 
use of a numerical grid adapted to the geometry of the 
system. We showed that this allows us to build high- 
accuracy initial data while keeping the number of grid 
points relatively low. 

Using the extended conformal thin sandwich formal- 
ism and fixing the initial state of the system through 
quasiequilibrium conditions, we obtained initial data 
whose precision is limited only by the small deviations 
from an exact equilibrium. As an example, we showed 
convergence tests for a sequence of equal-mass, irrota- 
tional BH-NS binaries, verifying the exponential conver- 
gence of our solver. Corotational and unequal-mass sys- 
tems lead to similar results. 

We also showed that with such accuracy we can resolve 
deviations from the point mass, circular orbit 3PN pre- 
dictions, and observe the influence of tidal distortion and 
eccentricity. 

Abandoning the assumption of conformal flatness, we 
generalized the method to construct binaries with a spin- 
ning black hole. Previously, initial data with a Kerr- 



Schild conformal metric was shown to be significantly 
inferior to conformally flat configurations [20|. Here, we 
showed that using a Kerr-Schild metric cut off at large 
distances from the BH allows reasonable precision to be 
reached — as in the case of BBH [2^ . We verified that 
with such a conformal metric we could construct a binary 
whose BH has a spin Jbh = —0.5 perpendicular to the 
orbital plane. 

Finally, we adapted a method designed for BBH [SOJ . 
and demonstrated our ability to significantly decrease the 
eccentricity of the binary initial data. 
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